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ABSTRACT 


Various  computational  aspects  have  been  investigated  to  numerically 
solve  charge-conservation  equations  and  Poisson's  equation  for  the  elec¬ 
tric  field  yielding  ion  and  electron  distributions.  These  equations 
were  derived  and  presented  by  the  Harry  Diamond  Laboratories  as  part  I 
of  this  study.  The  computational  aspects,  reported  herein  as  parr  II 
of  the  study,  include:  (1)  transformations  to  reduce  the  steep  slopes 
of  the  input  functions  and  to  simplify  the  solutions  of  the  equations; 
(2)  linearization  of  the  equations  to  permit  use  of  matrix  methods  in 
their  solution;  (3)  derivation  of  small-value,  asymptotic  solutions  to 
provide  starting  conditions  in  the  matrix  solution;  (4)  a  computer  pro¬ 
gram  listing,  description,  and  sample  output;  and  (5)  descriptions  of 
an  independent  check  solution  and  other  checks  to  confirm  validity  of 
the  results.  The  computer  program  is  written  to  accommodate  any  con¬ 
sistent  set  of  boundary  conditions.  Although  the  equations  are  linear¬ 
ized,  the  nonlinear  terms  are  approximated  in  a  way  to  insure  rapid 
convergence  of  solutions  to  the  exact  equations. 
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A  study  has  been  conducted  by  the  Harry  Diamond  Laboratories  of  ion- 
electron  distributions  for  a  hypersonic  vehicle  with  a  10-degree  semivertex, 
sharply  pointed  cone  at  Mach  8  and  10  for  sea-levi  1  flight.  Detailed  cal¬ 
culations  are  reported  by  Pollin1  as  pert  I  of  this  study.  This  report, 
which  represents  \  art  II~of-the  study-,  is  concerned  specifically  with 
computational  aspects  of  the  ion-electron-distribution  equations  -detailed 
in  part  I.1  The  equations  to  be  solved  are  presented  herein  to  describe 
the  transformations  required  to  simplify  their  numerical  solution.  Readers 
who  are  particularly  interested  in  the  derivation  and  precise  symbol  defi¬ 
nitions  should  refer  to  part  I.1 


Basically,  a  parabolic  system  of  five  first-order,  partial  differential 
equations  describe  the  phenomenon  (listed  as  equations  18  and  20  in  part  I) . 
Use  of  the  compact  "±"  notation  to  condense  two  equations  into  one  is  shown 
below. 


irvixj^l  _  VJAiSj 
L  P(n)  Jn  P  c 


i(x)  + 


(h)6(x)i 

P(n)  i 


J  +  ±  + 

6(x)P(n> - +  6(x)v(x,y)a  +  K(n)s<Ji 

en  n 


D(n)  +  DT(x,n) 

.  iflil  en 

8.35  ene  P(n)  (  8) 


Equations  (\)  and  (2)  are  each  two  equations — read  first  with  upper  and 
then  wiai  lower  signs  whenever  two  signs  appear.  The  dependent  variables 
+  -  +  - 

are  J,  J,  s,  s,  and  <J>v  :  the  independent  variables  are  0  <  x  <  00  and 
0  <  y  <  6(x)  w. n  defined  by 


so  that  0  <  T)  <  1  is  used  as  a  normalized  independent  variable.  This  is 

convenient  also  because  many  of  the  input  functions  are  ultimately  functions 

of  n.  The  other  symbols  are  either  constants  (p  ,  c  ,  e,  n  ),  or  other 

e  e  e 

+  _  i  ±  ± 

known  functions  [P(n),  v(x,y) ,  W(n,s,s),  S(x),  u(n) ,  D(n).  DT(x,n),  K(n) ] ; 

some  of  the  constants  and  functions  change  with  geometry,  velocity,  and 
altitude  of  the  vehicle. 


1Pollin,  I.,  "Ion  and  Electron  Distributions  in  the  Boundary  Layer  of 
Hypersonic  Vehicles  for  Chemical  Nonequilibrium  Flow — Part  I:  Aero¬ 
dynamics  and  Numerical  Results,"  HDL-TR-1565,  August  1971. 


Preceding  page  blank 


Subscripts  x  and  n  indicate  partial  derivatives  with  respect  to  these 

variables.  Note  that  does  not  appear  except  with  subscript  n;  <j>  is 

considered  as  a  basic  dependent  variable.  The  voltage  $  is  calculated 

once  <fc  is  obtained, 

n 

The  known  functions  are  sometimes  derived  from  formulas  and  sometimes 
are  calculated  in  point-data  form.  We  have 


S(x)  *  bx0,8,  b  constant. 

(5) 

u(n)  *  u.  constant, 

6  0 

(6) 

w(n,s,s)  =  KFN0(n)  -  B(n)ss 

(7) 

where  KFNO  is  used  for  K,<N><0>  and  B(n)  is  used  for  K  [n  /P(n)]2.  Both 

t  re 

K.  and  are  found  as  functions  of  the  known  temperature  function 
T  =  T(n): 


.,  SxlO*1 1 exp (-32 .500 /T) 

Kf  “  T^ 

K  -  3x10_3/T1 *  5 
r 

Also, 

.  11.600* 

K(n)  Tc^r°(T) 

DT(x,n)  *  0.02u(56(x)n9^8 
+  — 

and  D'VT)  and  D(T)  are  related  by 
D(T)  =  234D(T) 


(3) 


(9) 

(10) 


(ID 


Finally,  v(::,y)  must  satisfy 

[u  — —  n  1  /  81  +  xfvtx^l  =  0 

[u6  P(n)n  Jx  +  XlP(rl)Jy  U 


(12) 


with  v(x,0)  =  0. 
here, 


A  lengthy  computation  shows  that  in  the  case  considered 


v(x,y)  *  u6  ff'i)  (13) 

where  f(r,)  is  given  by 

f  (n)  «  0.8  n9/8  -  P(n)  J"  dt  (14) 

+  System  (1)  to  (3)  is  also  subject  to  boundary  conditions;  if  s(x,n), 
J(x,n)  ,  and  <J>  (x,n)  are  solutions,  we  have 
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s(0,n) 

=  s(0,n)  “0 

(15) 

•f 

8 (x ,0) 

■  s(x,0)  *  0 

(16) 

s(x,l) 

“  s(x,l)  ■=  ^(x.l)  -  0 

(17) 

THE  J' 

TRANSFORMATION 

A  convenient  simplification  to  system  (1)  to  (3)  results  from  trans- 

+  + 

formation  of  the  current  variables  J  to  J*  defined  by 
+ 

6(x)  (18) 


+ 

J' 


r  + 


_j_  .  sy(x,y) 

'  n  /  \ 


en 


P(n) 


By  taking  a  derivitive  with  respect  to  n* 


(J*) 


1  * 

+  Mi 

■)+  i 

V 

" 

[en 

P 

u 

L»  el 

n 

6(x) 


(19) 


and  substituting  into  equations  (1)  and  (2) ,  we  obtain  a  simpler  system, 
hers  expressed  in  terms  of  the  more  primitive  functions 


(J')  =  u  KOOJ2  |  .  L600J2  [KFNO(n)  -  B(n)ss] 
n  6  x  P(n)  n  pc  1  ‘ 

e  e 

+  u  -LiL*).]2  1/8  ± 

U6  P(n)  n  s> 

+ 

+  +11,600  s$n  +  P(n)J' 

8_  “  + 


n 


nn 


D(T)  +  DT(x,n) 

i sii.  en  (+  _ 

8.85  cne  P(n)  K  J 


(20) 


(21) 


(22) 


This  transformation  has  two  purposes:  (1)  to  eliminate  the  need  for 
derivatives  of  f(n)/P(n)  in  the  solution,  and  (2)  to  eliminate  the  singu¬ 
larity  6(x)/x  at  x  =  0.  This  term  now  appears  as  6?/x  =  b2x0,6. 


3.  SMALL-VALUE.  ASYMPTOTIC  SOLUTIONS 

Because  of  the  fractional  powers  of  x,  a  Taylor  series  solution  dees 
not  exist  at  x  =  0.  In  an  attempt  to  find  small-value  (of  x)  solutions, 
we  may  eliminate  the  fractional  powers  in  x  by  the  transformation 

t  -  x0 .2  or  t5  =  x  (23) 
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with 


dx 

dt 


=  St4*  , 


±  dt 
St  dx 


and 


6  (x)  =  bt4 

Equations  (20)  to  (22)  transform  to 

(H  -  v2t3  wvrr  [Kn,0(,') 

e  e 

+ 

*  11,600  U  +  P(n)3* 


-  B(n)ss]  + 


+ 

s 


nn 


D(t)  +  DT(t5,n) 

iO14  b2t8  ,+ 

"  8.85  eRe  P(n)  (S  "  S) 


5P(n) 


(24) 


(25) 


uAb2t4n1/8  § 

S  -1  (26) 


(27) 

(28) 


A  power  series  solution  now  takes  the  form 

S'  =  [  S  (n)t1 

J  A  1 


i*0 


(29) 


s  -  l  s .  (n)t1 
i-0  1 


(30) 


*n  "  J0(Vi(n)tl 


(31) 


By  substituting  equations  (29)  to  (31)  into  equations  (26)  to  (28)  and 

equating  coefficients  of  powers  of  t,  ordinary  differential  equations  are 

found  for  the  coefficients.  [This  procedure  first  requires  a  power  series 

+ 

representation  for  1/ (D  +  D^,) . ]  The  smallest  set  of  nonvanishing  coeffi¬ 
cients  occurs  for  subscript  8  in  equations  (29)  and  (30)  and  subscript  16 
in  equation  (31).  We  have 


(S8)r 


P  C 
e  e 


KFNO(n) 


(32) 


(s8)  =  J8 

n  D(T) 


10 


14 


16*n  '  ~  8.85  ene  P(n) 


(s8 


-  s8) 


(33) 


(34) 


10 


XV 


subject  to  the  boundary  conditions  of  equations  (16)  and  (17);  that  is, 

sg(0)  *  sg(0)  =■  t0(i)  -  ss(l)  -  (4>ri)  16 *  0*  Thus,  we  observe  that  for 
small  t  corresponding  to  small  x, 


j’(x,n)  *  JeOrOx1*6 

(35) 

4  +  l  r 

s(x,n)  *>  S8(n)x1,a 

(36) 

4>n(x,n)  *  (^n^ is ^n^x3 *2 

(37) 

This  analysis  also  suggests  that  <j>n  is  essentially  decoupled  from  the 

4-  _  + 

equations  for  s  and  s  for  small  x.  The  nonlinear  term  $-5  is  negligible. 

In  addition,  the  "+"  and  equations  are  decoupled,  since  the  ’onlinear 

term  ss  is  negligible.  These  results  are  almost  independent  of  the  boundary 
conditions  end  must  be  considered  when  sets  of  boundary  conditions  are 
desired  other  than  those  given. 


4.  THE  z,q  TRANSFORMATION 


The  independent  variables  n  and  x  are  inconvenient  variables  t£  compute 
with.  At  x  =  0,  we  have  from  equation  (?6)  that  sx  *  0,  but  that  s^  and 
higher  derivatives  are  infinite.  The  finite^difference  scheme  used  to 
approximate  sx  from  a  sequence  of  values  of  s  is  subject  to  large  error 
when  sxxx  is  large.  Therefore,  it  is  advantageous  to  use  smaller  Inter¬ 
vals  near  x  *  0  than  are  otherwise  used.  This  situation  is  reinforced 
by  the  expected  behavior  of  sx  +  0  as  x  +  the  intervals  in  x  can  then 

be  taken  further  apart.  To  accommodate  both  situations,  we  compute  in  a 
new  variable  z  defined  by 


z  «  x0'8  or  x  =  z1*25  (38) 

ipostly  to  simplify  checking  small-value  results  (s  ~  zO  and  provide  for 
szz  to  be  finite.  The  equations  are  changed  by  replacing  x  by  z3  *25  and 


±  ±  dz  0-8sz 
Sx  Sz  dx  z0*25 


(39) 


A  similar,  but  more  severe  situation  exists  for  the  n  independent 
variable.  The  driving  function  KFNO(n)  is  monotonically  decreasing  and 
is  almost  an  impulse  function.  For  the  Mach  8  data,  KFNO  diminishes  more 
than  three  orders  of  magnitude  in  the  range  0  <  n  <  0.01.  It  is  essential 
in  any  finite  difference  scheme  approximating  derivatives  with  respect  to 
n  to  have  many  intervals  for  small  r\.  A  transformation  that  automatically 
increases  the  number  of  r)  intervals  for  small  n  is  given  by 

q  =  in( n  +  a)  (40) 
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where  a  >  0  is  chosen  to  satisfy  any  particular  requirement.  In  this  case, 
a  was  chosen  so  that  in( a  +  0,0001)  -  £n( a)  was  equal  to  1/400M  of  the 
total  range  of  q.  Equivalently,  this  assured  that  when  400  equa.'  intervals 
were  taken  in  the  q  direction,  the  first  value  of  q  corresponded  to  n  *  0, 
and  the  second  to  n  =  0.0001.  The  procedure  was  motivated  by  the  fact  that 
many  of  the  input  functions,  including  KFNC(q) ,  were  constant  for 
0  £  n  <  0.0001  and  dropped  of  steeply  thereafter.  A  special  computer 
program  computed  6t(a)  =  -4.7939598. 


From  equation  (40) ,  we  have 


q 

n  *  e  -  a 

ia  = 

dq 


(41) 

(42) 


This  transformation  is  made  by  replacing  n  by  e1*  -  a  [eq  (41) ]  and  replacing 
any  variable  rq  by 


In  practice,  the  n  notation  was  retained,  except  for  derivatives  with 
respect  to  q.  All  functions  were  computed  as  functions  of  n  as  found  from 
equation  (41)  for  any  q.  The  transformed  equations  [from  (20)  to  (22)] 
now  appear  as 


m  u  i 

eq  el  z>'2S  P(n)  Sq 


x2/1.25\  , 

— - L  [KFN0(n)  -  B(n)ss] 


+  0.8u,. 
o 


62(z1.25)nl/8  ± 
Z0.2S  p(n)  Sz 


(44) 


sq  +11,600  -4-4  S(j-n  +  P(n)j' 

"Tf  =  - + - - 

D(T)  +  DT(z1-25,n) 

^n^q  _  _  IQ1'*  62 (z1  *25)  .+ 

e*!  '  "  8.85  ene  P(n)  (S 

with  6(z1 *25)  *  bz. 


s) 


(45) 

(46) 
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5.  NUMERICAL  SOLUTION  METHOD — LINEARIZATION 


A  "marching"  technique  is  used  to  solve  system  (AA)  to  (A6)  numerically. 

+  + 

Given  solutions  at  z  *  zj-i»  i*e.,  I(zj_1,q),  J'(zj_i»q),  and  4>n(zj_],q), 

ordinary  differential  equations  in  q  are  found  at  z  «  z*  *  zi.i  +  Az  by  an 

±  J  J  ± 

appropriate  finite  difference  approximation  for  §z  in  terms  of  s(zj,q), 
s ,q) ,  etc.  To  simplify  notation  hereafter,  we  refer  to  a  function  at 
fixed  Zj  with  subscript  j;  i.e.,  s(zj,q)  e  I j  ,  understanding  that  it  is 
also  a  function  of  n. 


The  marching  technique  requires  the  solution  of  a  two-point  boundary- 
value  problem  with  nonlinear  differential  equations.  Although  these  can 
be  solved  directly  by  iterative  processes,  faster  matrix  algorithms  can 
be  used  for  linearized  equations,  which,  at  the  same  time,  do  not  materi¬ 
ally  affect  convergence  characteristics.  Accordingly,  system  (AA)  to 
(A6)  is  linearized  by  deriving  differential  equations  for  changes  in 
i  — 

Sj.j  ,  Jj-1»  and  (♦  Defining 


■  gJ-l + 

(A7) 

+  +  + 

■  5j-I  +  "j 

(A8) 

"  (Vj-l 

+ 

(A9) 

note  that  system 

(AA)  to  (A6)  holds  at  both  values  of  z^_i  and  Zy 

We 

substitute  equations  (A7)  to  (A9)  into  equations  (AA)  to  (A6)  and  subtract 


from  the  resulting  equations  a  similar  set  of  equations  obtained  when 

X  ± 

s  ,  J'  ,  and  (<j>  ).  ,  are  substituted  into  equations  (AA)  to  (A6) .  This 

J—  *  J —  1  h  J  — *  4  +  4.  ± 

yields  differential  equations  for  As^ ,  AJ^,  and  (A<t>n)j  with  s^  i  *  i> 
and  (<J>ri)j_1  as  additional  input  functions.  The  final  equations,  however, 

depend  upon  how  th  derivative  i  is  treated,  as  well  as  the  nonlinear 

±  +  -  z 
terms  As.(A<|>  ),  and  As. As.. 

J  Vj  j  j 


A  backward  difference  approximation2  is  used  for 
at  z  =  z, 


r 


± 

8 


3*j  - 


+  + 
Asj_j  +  s 


2Az 


“ —  +  0(Az2) 


+ 

sz  when  substituting 


2Kcpal,  Z.,  "Numerical  Analysis,"  John  Wiley  &  Sons,  New  York,  1955, 
pp.  515-518. 
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s 

z 


3As j  —  As j_j 
2Az 


+  0(Az2) 


(50) 


+  +  + 

where  s  =  s  -  s.  ,  a  quantity  that  must  be  stored  during  the 

J  —  1  J*"^ 

computation.  The  appendage  0(Az2)  indicates  the  approximation  is  of 
second  order  in  Az. 


Equation+(50)  is  not  valid  when  z  *  Az,  i.e.,  for  the  first  yalue  of 
z,  because  A§0  is  nonexistent.  Here,  ve  use  the  condition  that  §2  *  0  at 
z  =  0.  We  may  then  show  that  the  approximation  becomes 


2As  i 


Az 


+  0(Az2) 


(51) 


A  central  difference  approximation  is  used  for  s  when  substituting 
at  z  *  2j-i 


+  + 

J"^Azlr~  +  °^zl) 


±  ± 
As^  +  As 


2Az 


-  +  0( Az2) 


(52) 


When  z  =  Az,  we  tcke  =  0  at  z  s  0,  since  it  is  a  boundary  condition. 
Similarly,  it  is  shown  in  appendix  A  that 

A®j(AVj  =  °  +  °(Az2)  (53) 

ASjAs ^  =  0  +  0(Az2)  (54) 

The  substitution  of  equations  (50)  to  (54)  yields  second-order  computations 

of  As,,  AJ , ,  and  (A<p  )  which,  in  turn,  produce  linear  convergence  in  z  for 
+  ±  J  3  r'J  +± 

s,  J',  and  .  (If  Az  is  halved,  each  As,  AJ ,  and  A^  is  produced  with  one 
fourth  the  error,  but  there  are  then  twice  as  many  intervals  to  sum  to  reach 
a  given  z.  Hence,  if  z  is  halved,  the  final  error  is  only  halved.)  How¬ 
ever,  it  was  found  that  the  error  introduced  by  equations  (53)  and  (54)  was 
greater  than  that  of  equations  (50)  to  (52).  Accordingly,  equations  (53) 
and  (54)  were  changed  to  (see  appendix  A) 


As .  (A<p  ).  =  As,  (A  4>  ).  +  0(  Az3) 

3  n  j  j-r  n  j 

(55) 

As. As.  =  As.  , (As) .  +  0( Az3) 

3  3  3”1  3 

(56) 

=  As^_j  (As) j  +  0( Az3) 

which  still  are  linear  terms  in  subscript  j  variables. 


14 


The  final  system  of  linear  equations  that  must  be  solved  for  each  j  is 
given  by  (the  subscript  j  is  hereafter  omitted  to  emphasize  the  dependent 
variables;  As  £  As^ ,  etc.): 


,-q  d£A|Q. 


dq  D(n)  +  D„(zJ-25fi 
*  J 


{  +K(n)[(lJ_1  +  AlJ_1)AcJ»n  +  (^Jj.jAi 

A¥j-i(V3-i  i 

D(n)  +  DT(z^25,n) 


f  P(h)[aj  -  (J‘)  T 


^  1  D(n)  +  DT(z^25,n) 


&]  Aic 


e-q  <L£Ml  _  iinl  r  +  L  _  KfflQtq)^2 

6  dq  U6  P(n)  [7™*  eq  +  ^pT2Tj  eq  J 


Az  P(n) 


-q^KJ_  lO^f^r  2/  + 


Fwhz(‘" '  +s{2<Vi '  Vi>) 


where 


AD„  =  0.02uxAzn 
T  o 


A  ^  j  J”* 

A 1-1.251  =  _1.25  "  _1.25 

1  '  Zj  Vl 

A62  =  S2.  -  62 

j  j-1 


If  As0  is  initially  set  egual  to  0,  the  correct  formulas  are  obtained  when 
j  =  1  if  the  term  362/z°,<!-5  is  changed  to  A62/Zj’25  in  equation  (58). 


It  should  be  emphasized  that  quadratic  convergence  In  Az  can  be  ob- 

±  - 

talned  If  all  approximations  are  such  that  the  As,  AJ,  and  A$n  are  com¬ 
puted  with  order  Az3.  Although  this  was  not  done,  appendix  B  indicates 
the  changes  required. 


5.1  Matrix  Formulation 


If  a  grid  of  N  +  1  equally-spaced  points  is  placed  along  q,  starting 
at  q'j  corresponding  to  n  *  0  and  ending  at  q^+1  corresponding  to  n  ■  1, 

we  may  correspond  to  these  points  5(N  +  1)  functional  values  we  seek. 

Thus,  for  every  j,  we  seek  for  all  c. . ,  1  <  i  <  N  +  1  the  values  As.,  ,  As. , 
±  1  11* 

AJ^ ,  and  (A$  )  .  Approximate  values  may  be  found  with  second-order  error 

in  Aq  by  satisfying  equations  (57)  to  (62)  at  the  half-interval  stations 
by  using  the  approximations 


(Vi+l/2 

ri+l/2  " 


—  +  0(Aq2)  1  <  i  <  N 

+  0(Aq2)  1  <  i  <  N 


(63) 

(64) 


where  r  is  any  required  variable  or  stored  variable  at  j  -  1.  This 
results  in  5N  linear  algebraic  equations,  which,  coaled  with  the  five 
boundary  conditions  of  equations  (16)  and  (17),  nay  be  uolved  for  the+ 
5N  +  5  functional  values.  By  incrementing  z j ,  storing  the  required  As 
at  each  i  (to  be  the  new  set  of  Asj_j),  and  mechanizing  equations  (47) 
to  (49),  the  procedure  may  be  used  to  reach  any  desired  value  of  z. 


5 . 2  Computer  Solution 

The  computer  program  is  written  to  accommodate  any  set  of  feasible 
boundary  conditions,  since  this  was  one  uncertain  aspect  at  the  time  of 
development.  It  is  important  to  prearrange  the  equations  so  that  their 
final  matrix  form  has  a  matrix  coefficient  that  is  block-tridiagonal 
(each  block  a  5*5  square  matrix) ,  in  order  to  use  an  efficient  algorithm3 
*or  their  solution.  A  difficulty  in  doing  this  is  the  uncertainty  of 
where  the  boundary  conditions  will  be  imposed.  There  must  be  five  con¬ 
ditions,  but  m  of  these  will  be  at  n  =  0  and  5  -  m  at  n  =  1. 


3Isaacson,  E.  and  Keller,  H.B.,  "Analysis  of  Numerical  Methods,"  John 
Wiley  &  Sons,  New  York,  1966,  pp.  58-61. 
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The  method  uatd  is  described  with  the  aid  of  figure  1,  where  the 
scheme  is  drawn  fcr  the  boundary  conditions  expressed  by  equations  (16) 
and  (17).  An  array  IDX  is  read  in  with  values  0  or  1.  As  shown  in 
figure  1,  IDX(l)  corresponds  to  s,  etc.  If  IDX  *  0,  then  the  finite 
difference  equation  obtained  at  station  i  *■  1/2  (the  vertical  arrow) 
is  used  in  the  i-th  group  of  five  equations  (the  diagonal  arrow) .  If 
IDX  -  1,  the  finite  difference  equation  obtained  at  station  i  -  1/2 
is  used  in  the  i-th  group  of  five  equations.  In  the  first  group  of 

■f  - 

five  equations,  there  are  no  AJ  and  AJ  finite  difference  equations; 
these  are  replaced  by  the  boundary  conditions  (in  this  case.  As  *  0  and 
As  «'0).  In  the  last  group  of  five  equations,  there  are  no  a£,  oi,  and 
A^  finite  difference  equations;  these  are  also  replaced  by  the  boundary 
conditions  As  *  0,  As  *  0,  and  A^  ■  0.  Note  in  this  case  that  1DX(1)  *  1 
and  IDX(3)  *  0  could  have  been  used  because  of  the  symmetry  in  the  bound¬ 
ary  conditions;  similarly,  IDX(2)  *  1  and  IDX(4)  *  0  could  have  been  used, 
"he  final  coefficient  matrix  appears  as 


where  each  nonzero  entry  is  a  5><5  matrix.  If  IDX(L)  =  0,  nonzero  elements 
are  contributed  to  Aj^  and  cn  line  L.  If  IDX(L)  -  1,  nonzero  elements 

are  contributed  to  and  on  line  L.  The  procedure  to  solve  the 

resulting  equations3  requires  the  storage  of  N  5x5  T  matrices  and  5(N  +  1) 
y  elements.  Thus,  the  limit  of  N  is  about  400  on  the  HDL  7094  without 
external  storage.  A  more  complete  description  and  listing  of  the  code 
is  given  in  appendix  C. 

3Isaacson,  E.  and  Keller,  H.B.,  "Analysis  of  Numerical  Methods," 

John  Uiley  &  Sons,  New  York,  1966,  pp.  58-61. 
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Figure  1.  Scheme  for  matrix  generation. 


5 . 3  Solution  Checks 


A  separate  program,  COMPF  (not  included  here),  computed  f(n)  using 
equation  (14)  from  the  P(n)  data  fed  in,  and  punched  the  data  out  in  the 
proper  format.  It  utilized  a  standard  integration  routine;  the  final 
results  were  checked  approximately  by  Simpson's  rule  for  the  Mach  8  data. 
Thereafter,  the  program  was  assumed  correct  for  other  data. 

Another  separate  program,  FTEST  (also  not  listed  here),  checked  out 
the  FUNCT  subroutine  listed  in  appendix  C.  This  was  important  to  insure 
that  all  functions  of  n  were  properly  generated  and  stored  for  use  in  the 
main  body  of  the  program.  Included  was  a  printout  utilizing  the  stored 
functions.  Points  were  spot  checked  in  each  column  .o  insure  reasonable¬ 
ness.  The  printout  for  the  Mach  10  data  is  shown  in  table  I. 

In  addition  to  the  printout  of  all  functions,  graphs  of  the  functions 
were  plotted  against  n  using  HDL's  CalComp  plotting  equipment  by  another 
specially  written  subprogram  GRAPHS.  In  this  way,  gross  errors  or  any 
disturbance  from  the  quadratic  interpolation  of  the  input  data  could  be 
easily  detected.  Examples  of  the  graphs  are  shown  in  figure  2. 

The  major  check  on  the  validity  of  the  results  was  a  separate  inde¬ 
pendent  program,  POLCHK,  which  solved  system  (57)  to  (59)  directly  with 
the  aid  of  a  d^  irferential-equation-solving  subroutine  FN0L2.  Solutions 
of  two-point  boundary-value  problems  require  no  iteration  with  linear 
differential  equations.  Let  the  vector-matrix  differential  equation  be 

A(n)  ^  =  g(n)  (66) 

+  -  +  — 

with  xi  =  As,  x2  =  As,  X3  =  AJ,  X4  =  AJ,  and  X5  =  AS^.  We  start  at 
n  =  1  and  run  solutions  with  decreasing  n,  using  the  known  initial  (at 
q  =  1)  conditions  xj(l)  =  0,  X2(l)  =  0,  and  Xs(l)  =  0.  The  other  con¬ 
ditions  vary  in  each  of  the  three  separate  solutions: 


A(n)  =  g(n) , 

P3  =  df ,  P4  =  0 

(67) 

A(n)  =  g(n) , 

q3  =  0,  q4  =  d2 

(68) 

A(n)  ~  =  g(n) , 

r3  =  dl >  r4  =  d2 

(69) 

where  pR(l)  =  qk (1> 

=  r^(l)  =  x^O)  =  0  when  k  =  1,  2,  or  5. 

Because  of  the  superposition  theorem,  it  is  easily  shown  that 


x  *  ap  +  bq  +  cr,  (70) 

if  and  only  if 
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a  +  b  +  c  =  I,  (71) 

and  the  boundary  conditions  are  satisfied: 

As (0)  =  xj(0)  *  api(0)  +  bqi(0)  +  crj(0)  =  0  (72) 

As (C)  =  X£ (0)  «  ap2(0)  +  bq2(0)  +  cr2(0)  *  0.  (73) 


Equations  (71)  to  (73)  can  be  solved  for  a,  b,  ana  c.  The  initial  con¬ 
ditions  and  d2  can  be  arbitrarily  nonzero,  but  in  practice  they  are 
chosen  so  that  r(n)  is  close  to  x(n);  i.e.,  c  will  be  very  nearly  equal 
to  1.  This  minimizes  roundoff  errors  since  the  magnitudes  of  a,  b,  and 
c  can  otherwise  be  quite  large  but  of  opposite  sign.  The  program  was 
written  so  that  di  and  d2  were  iterated  on  successive  runs  h  and  h  +  1 
before  z  was  stepped  up: 

(di)h+1  *  (di)h(a  +  c)  (74) 

(d2)h+1  =  (d2)h(b  +  c)  (75) 


The  iteration  was  stopped  when  the  magnitudes  a  and  b  were  appropriately 
small  enough.  Theoretically,  c  at  step  3  should  equal  c  at  step  2  exactly, 
but  because  of  roundoff  errors  and  the  inability  to  exactly  solve  differ¬ 
ential  equations  numerically,  the  iteration  index  occasionally  went  to  4. 

A  fourth-order  Runge-Kutta  integration  scheme  was  used  to  compute  the 
values  of  variables  at  fixed  points.  These  points  were  sufficiently  close, 
so  that  error  in  the  n  direction  was  essentially  negligible.  This  check 
showed  that  the  finite  difference  scheme  of  section  5.2  was  working  cor¬ 
rectly,  and  that  the  error  introduced  in  the  n  direction  was  less  than 
2  percent  when  using  400  intervals  (for  Mach  8  data) . 


By  obtaining  runs  at  a  given  Az,  Az/2,  and  Az/4,  the  linear  nature 
of  the  convergence-  in  the  z  direction  was  verified.  In  addition,  the 
error  using 


Az 


100°  *  P  _  zmax 
50  =  50 


was  less  than  5  percent  (for  Mach  10  data) .  Such  a  run  corresponded  to 
about  8  minutes  of  IBM  7094  CPU  time. 
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Appendix  A.  Order  of  Errors  In  Nonlinear  Approximations 
Since 


+  1 

±  ±  ,  9s  ,  ,  ,  1 

6‘Vl+  ^L<2'Vl)+2 


we  have 


'  Jj-i  ’  slj  ‘ 


^  ,  .  i 

-  j->4z  +  5 


0 

32s 

3z2 


Az*-  +  . . 


/j~l 


Similarly, 

(AVj 


J-l 


Az  +  | 


a2<J>r 


s72- 


Az2  + 


J-l 


so  that 


4V4Vj 


+ 

3s 

l**n 

3z 

j-i 82 

Azz  +  • . • 


j-1 


=  0  +  t>(Az2) 


H-i 


9s 

3z 


+ 


•  •  • 


*  0  +  0(Az2) 


This  implies  that  ignoring  the  nonlinear  terms  simply  introduces  a  second- 
order  error.  On  the  other  hand,  we  also  have 


+ 

Ai 


J-l 


+ 

3s 

3z 


j-2 


Az  +  j 


32s  1 


Az21  + 


so  that  by  subtracting 


This  leads  to 

As^  =  As^^  +  (?(Az2) 
and 

A$j(A*n)j  =  As  (A*n)  +  0(Az3) 

Equation  (56)  follows  similarly.  Since  ASj_j  must  be  stored  for  the 
derivative  approximation  equation  (50)  anyway,  no  additional  storage 
is  required. 
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Appendix  B.  Formulas  for  Quadratic  Convergence  in  z 

+  ± 

The  first  set  of  As,  AJ,  and  A<jn  can  be  computed  with  a  second-order 
error,  but  each  set  thereafter  should  be  computed  with  third-order  to 
obtain  quadratic  convergence.  +No  change,  therefore,  occurs  for  j  =  1; 
equation  (51)  applies  for  the  s2  approximation  at  z  =  Az. 

For  j  =  2,  i.e. ,  z  =  2Az,  a  special  computation  must  be  made,  taking 
into  account  the  stored  AS]  and  the  fact  that  sz  =  0  at  z  =  0.  Note  that 
we  need  to  estimate  sz  at  z  =  Az  as  well  as  z  =  2Az.  If  z  =  Z](=  Az) , 
we  set  up  a  Taylor  expansion  of  s(z)  about  z1# 

S  =  S]  +  S  J (z  -  Z])  +  J  s"(z  -  Z])2  +  ~  s'"  (z  -  z  x ) 3 


zl)4 


+  .  •  . 


+  S  j  (z  -  Zj)  +  -J  s'"  (z  -  ZX)2  +  s\"(z  -  Zj): 


+  •  •  • 


At  z  *  Zj  -  Az  (=  Zq  =  0) ,  s  *  0  and  s'  *  0.  Also,  at  z  ~  Zj  +  Az 
a  *  s2.  Thus, 

0  *=  s0  =  sj  -  s{Az  +  y  s"Az2  -  j-  s"'Az3  +  s'1",Azlt  +  ... 

0  =  sj  -  s'j'Az  +  ^  s'J'Az2  -  j  s""Az3  +  .  . . 

s2  -  sx  +  s}Az  +  j  s"Az2  +  ^  s"'  Az3  +  s'^Az4  +  ... 


=  z 


2 » 


Eliminating  s"  and  s'"  from  the  equations  yields 
s'  =  Asl  +S2.  -  IT  a,3 


4Az 


tzt 


Az3  + 


or,  in  terms  of  As2  and  Asj, 

5As,  +  As,  . 

si  — te — 1  +  0(Az  > 

since  s2  -  Sj  =  As2  and  Sj  =  AS]. 

To  find  s2  in  terms  of  Asj  and  As2,  a  Taylor  series  about  z 
constructed: 


z2  is 


s  =  s2  +  s2(z  -  z2)  +  j  s'2(z  -  z2)2  +  — 


s  2  '  o  2 

i 


s2"  (z  -  Z2) 


+  -S'2"\z-z2)4  + 
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with 


s'  =  sj  +  s  2  (z  -  z2)  +  y  s!,”  (z  -  z2)2  +  ~  s’2”'(z  -  z2)3  +  ... 

At  z  =  0  (z  -  z2  =  -  2Az)  ,  s  =  s'  =  0,  and  at  z  =  Zj,  s  =  Sj.  Thus, 
0  =  s2  -  2s2Az  +  2s2Az2  -  |  s2'  Az3  +  st>”’Az4  +  ... 

0  ^  s2  -  2s'2Az  +  2  s2"  Az2  -  ^  s2'"Az3  +  ... 

S1  53  s2  ”  s2Az  +  2'S2Az2  -  ~  s2"Az3  +  —  s2'"Az4  +  ... 

Eliminating  s2  and  s2'  from  the  equations  yields 

2s-  -  'vs,  i  , 

s'0 - V- - L  +  £S’2"'Az3  +  ... 


Az 


In  terms  of  As2  and  Asi, 


2As,  -  2As,  ,  . 

*  - \-z - L+0( Az3) 


A  projected  ASj  for  use  in  the  nonlinear  approximations  (55)  and  (56) 
can  also  be  found  by  eliminating  s2  and  s2'  from  the  above  equations. 
We  have 

2 

s2  =  As  i  +  —  s'”  Az3  +  . . . 


or 


As2  =  3Asj  +  0(Az3) 

From  appendix  A, 

As2(A4>  )2  =  SAi^A^n^  +  ^(Az4) 


>2'-.n- 

As2As2  =  3AsiAs2  +  0(Az4) 
.+ 


or 


=  3As2Asj  4-  0(Az4) 


For  j  >  3,  we  may  use  standard  formulas2 


M-i 


S1~3  -  6s.1-2  *  3sj-l  *  2s1 
6Az 


+  0(Az3) 


2Kopal,  Z.,  "Numerical  Analysis,"  John  Wiley  &  Sons,  New  York,  1955, 
pp.  515-516. 
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Appendix  C.  Program  Description  and  Listing 


in  the  listing  of  the  computer  program  (presented  on  pages  30-40) ,  the 
main  program,  P0LMD2,  utilizes  the  following  subroutines  to  simplify  the 
flow. 

FUNCT — Reads  in  the  main  data,  computes  some  needed  constants  and  all 
needed  functions  of  n  at  half-integer  stations  by  quadratic  interpolation. 

TERP2N  (not  listed) — Called  by  FUNCT  to  perform  quadratic  interpolation. 

INITLZ — Initializes  all  variables  and  computes  other  needed  constants. 

The  subroutines  FUNCT,  TERP2N ,  and  INITLZ  were  overlayed  by  subsequent 
subroutines ,  since  they  ore  needed  only  at  the  start  of  the  run. 

XSTEP — Steps  z  and  computes  all  constants  that  are  functions  of  z  only. 

EMPTY — Initializes  the  A^,  B^,  and  q  matrices,  and  the  q  vector  to  C. 
Note  that  the  A^  matrix  is  called  X  in  the  program,  the  Bj  is  called  BB, 
and  the  q  matrix  is  adjoined  with  the  q  vector  to  form  a  5*6  matrix 
called  D.  (This  permits  the  later  computation  of  A^C  and  Aj*q  in  one 
operation. ) 

EMPTY  starts  the  minor  i  loop,  the  purpose  of  which  is  to  eliminate 
the  i-th  block  of  five  equations 


.+ 

+ 

/  + 

ASi 

\ 

ASj 

ASj  . 

As  j  \ 

J  i 

Asj 

I  • 

+ 

+ 

1  + 

4Jj  | 

i  +  A,  AJ, 

1  j 

+  q  j  AJj  -  q 

“j  j 

AJj 

1^  Uj 

,-i  i 

\  (A<Mj  i+1 

(Bj  =  0  for  i  =  1  and  q  =  0  for  j  =  N  +  1)  where  the  subscript  j  indicates 
that  z  =  jAz. 

FIRST,  SECOND,  THIRD,  FOURTH,  FIFTH— are  entries  into  EMPTY  that  fill  the 
X,  BB,  and  D  matrices  line  by  line.  Note  that  single  subscript  arithmetic 
is  used  to  conserve  computer  time  and  core. 

ENDO,  END1 — are  called  when  i  =  1  and  N  +  1  to  take  the  boundary  conditions 
into  account. 

BOUNDO,  B0UND1 — define  the  boundary  conditions  at  n  =  0  and  r,  =  1.  This  is 
the  only  subroutine  the  user  must  write  if  his  boundary  conditions  differ 
from  equations  (16)  and  (17).  [Actually,  the  supplied  B0UND1  cor tained  the 
condition  that  <tiri(x,l)  =  P^;  Pj  was  read  in  to  be  0  to  satisiy  equation 
(17).] 
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XDCORR — corrects  the  and  arrays  by  the  matrix  multiplication: 

“  A^  -  ,  i  *  2,3, .. . ,N+1 

fj,  •  fj  ”  ^i^i-l  »  ^  “  2,3,...,N+1 
as  required  by  the  algorithm.3 

To  speed  the  process,  only  the  possible  nonzero  elements  of  were 
used  and  every  term  was  written  out. 

GELG  (not  listed)  —  computes  A^D  by  Gaussian  elimination  using  complete 
pivoting. 

The  first  five  columns  of  Aj*D  are  stored  as  the  array  (in  the  com¬ 
puter  Z  array)  and  the  last  column  of  AJ*D  is  stored  as  the  y^  array  (in 
the  computer  XX  array) . 

REST — does  the  rest  of  the  computation  in  the  i  loop: 

1.  Back  multiplies  to  find  the  **ariables. 

2.  Updates  all  needed  variable  functions  of  n  and  some  functions 
of  x. 

3.  Determines  if  a  printout  is  required  for  the  value  of  z  just 
completed. 

4.  If  a  printout-  is  needed,  computes  by  rectangular  formulas, 

±  ±  J 

retransforms  to  J  from  J’,  and  computes  other  functions  as 
needed  for  the  printout.  A  sample  printout  is  shown  in  table 
C-I  on  page  41. 

In  the  main  program,  the  unlabeled  COMMON  is  placed  in  equivalence 
with  the  Z  array  to  save  core.  All  input  data  and  some  variables  needed 
in  the  FUNCT,  INITLZ,  and  REST  routines  are  not  needed  during  the  i  loop 
and  can  be  shared  with  the  Z  array.  Liberal  use  of  labeled  COMMONS  served 
to  otherwise  transfer  data  from  subroutine  to  subroutine. 


3Isaacson,  E.  and  Keller,  K.B., 
Wiley  &  Sons,  New  York,  1966, 


"Analysis  of  Numerical  Methods,"  John 
pp.  58-61. 
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P0LM02 . EFN  -  SOURCE  STATEMENT  .  .IFN(S)_.--. - 


DIMENSION  Z(  10300 1  -  -  -•  -  —  - 

COMMON  ETAPTI 130) .PPTI 300),ETALG( 100 ) ,BOB( 200) ,TPT2 ( 100) , DPT ( 100) 

)  0(?>  ,F1(3>  ,PMJ(<*01  ) ,NS2,XN,XZ,XEND,XPR!NT,SNE,DY,Y0,ETAI,X:.Y. 

2  EXPO.uPLS.Al ,42, A3 , OELT AX , JOP , JOM, VSM, PHI Y, PSYSP .USYSM , SNM, 
s  TSYSP.TSYSM 

C'jmj,j/CQNSTI/XJ»OEUX,ZJ,0!25,OELJ 
COMMON /CONST 2 /PI ,  P2  »  P3, P4, P5  *PS  .  - 

COyMON/CONST3/  0CLZ,CU'L2,RH0ECE, V0X,ENEE.XZ12,XJ12,DELJ12,DELJX2 
1  , C 1  8  ) ,  V  ,N£ ,  CUEL »CDELt,ODCLTA, DELTA 

CG*N.UN/FUNl/OPI.US1400),OMlNUS(400>  ,EPKOO)  ,ETA98  1  400)  ,T29(4D0), 

-  1  214.0)  ,FPEY(43C),KfNn(4.?0),DY-Y(4C3>,Pi400),T29M(430> 

CDwv.uN/ FUN2/  B3,81,Krr.0D,KF  101 

COMMON/ FUN2/  ETAC(51),FO(al),OTPLO(51),FPOC51),PO(51)fOPLUO(51l#  - 
1  OPKlUi^l ),£TA9C(51> 

C OPrlO.N/  INTLZ/PHI JI401)  t  SPJI401  >  i  SMJ 1 401 ) « JPJl  401  >  *  JMJ  (401 )  • 

.  1  USMJ<401).DSMJ<4Gl),SPJH(4G3),SMJitI400),PHIJH(40C)  . . 

C0MK0N/MATI/X{25),U(3O),8bC’5) 

CQ-W>'0N/HAT2/  XXI  2005 ) 

.....  COMhOJ/lNOl/  I.  J.M,  II,JJ,K  - 

COMMON/ 1 NOX/  1DX16) 

COM. MON/ NUMB/  N , N W , N P l , N C  »  NC  J UN T 

EQUIVALENCE  1Z1  1),  ETAPTI  »)*  _ _  ... 

HEAL  JOP. jnM,xFNG,JPJ,JMJ,NE,NE2,KFN00,KFN0l 
C  COMPUTE  FUNCTIONS 

CALL  FUNCT  -  ...  - -  ... 

C  INITIALIZE  ALL  VARIABLES 
CALL  INITLZ 

C  ROUTINES  FUNCT,  INITLZ,  AND  TEPP2N  ARE  NO  LONGER  NEEDED.  . . . 

C  MAJOR  Z  LOOP 

00  9V9  J  «  1,M 

c-  initialize  x  functions.  ...  — - - : - 

CALL  XSTEP 
C  MINOR  cTA  loop 

00  900  I  «  l.NPl  .  - - - - - - 

C  EMPTY  A no  FILL  X,  0,  ANO  08  MhTRICES  WITH  A,  C  AND  F,  AND  B. 

CALL  EMPTY 

.  IFll.tU.l)  GC  TO  10  - - - - - 

I  F ( I  .C0.NP1 )  GO  TO  20 
CALL  FIRST 

CALL  SECOND  . .  .  _  _ 

CALL  THIRD 
CALL  FOURTH 

_ CALL  FIFTH - 

GO  TO  i 
10  CALL  ENDO 

20  CALL  EN01 

C  CORRECT  FOR  8  WHEN  I  NOT  EQUAL  TO  1 

3  CALL  XOCORRIZ ) .  —  - - 

C  ^EAOY  rCf  INVERT 
150  CALL  GLLG1 0, X ,5, 6 ) 

C  TRANSFER  GAMMA  MATRICES  ANO  Y  VECTOR  TO. STORAGE.. _ _ 

IF  I  I .tO.NPI )  GO  TO  200 
JJ  =  2  5  *  1 1 - 1 ) 

O'J  iSJ  II  =  1,25  . .  .  - 


NOT  REPRODUCIBLE 


A  *  1 1 

♦  JJ 

l&d 

l  IK)  « 

OMI) _ 

2C3 

JJ  *  5MI-1) 

00  210  I!  *  1,5 

. 

K  «  JJ 

♦  11  -  ... 

210 

XX (Kt  * 

0  ( 1 1  ♦  .’5 1 

C  END  OF  £T 4  LOOP. 

900  CONTINUE  —  .  _ — 

C  UPDATE  AilO  WRITE  IN  REST  ROUTINE 
CALL  RcST ( Z  I 

0.  ENO  OF  l  LOOP _ _ _  . 

‘>99  CONTINUE 
STOP 

..  -  ENC  - 


- FUNCT. _ --..EFN  SOURCE.  STATEMENT. _ - _ IFN(S»_.-_ 


SUBROUTINE  FUMT .  .  .  .  .....  _ _ 

COMMON  ETAPT ( 130) , PPT ( 300),ETALC( 100),B08(200),TPT2( 100) .DPTIIOO) . 

1  0(2),FU3),PHJ<*0l),NE2,XN,XZ,.'<END,XPRINT,ENE,DY,Y0»ETAI,XI,V, 

2  S<P0,0PLS,A1,A2,A3,UE1.TAX,JQP, J0M»V5M,PHIY,DSYSP.0SYSM,SNM, 

3  TSYSP.TSYSM 

COMMON/ C0NST3/  DELZ,CDEL?,RH0ECE,VDX,E.NEE,XZ12,XJ12,DELJ12,DEUX2, 

t  OZ232«C(8)«V ,NE  «CUEL«C DELE, OOEIT A, DELTA . .  .  . - 

COMMON/ FUN l / DPIUS ( AOO )  .UMINUSI 400 ) , EPJ 490 ) »  ETA98 1  400 )  »T29( 40G )  • 

•l  814 Ju)*FPSY( 490)#KFNQ 1400 )»DYEY(400)tP( 400 ) , T29MI 400) 

COMMON/ FUN2/  09, 31»KFN09,KFN0t  .  _ 

CUMM0N/FUN3/  ETAO( 5 1 ) » F0( 51) .OTPLOI 51 ) ,FPO( 51 ) .P0( 51 ) .DPLUO! 5 1 ) , 
l  0PMI0(51),ETA?0(5l) 

-  -  COMMON/  INDl/  I ,  J  .  M  .  lltJJ(K  ..  _ _ _ _ 

COMMON/ NUMB/  N, NW.NP  l, NC .NCOUNT 
P.EAL  NC.NE2,  JO0.  JOM.KFNO.KFNOO.KFNOl 

READ (5 • 1039)  V,CUEL,RH0ECE,N£,CDEL2  _ 2 _ 

10U0  FORMAT ( 6E13. 5 ) 

REALMS,  jO)  I  I  ,  J  J  t  K 

.  30  FORMAT (315)  . - 

REALMS,  1000)  (  ETAPT(l),  1-1,11) 

REAOIS, 1000) (  PPT11J, 1-1,1!) 

J  *  II  ♦  l  - - 

N  -  II  ♦  II 

REALMS  ,  1000 )  (  PPT(H,!»J,N) 

-  -  J  *  N  ♦  1  —  - - - 

N  -  N  ♦  II 

REAOIS. 1009) I  PPT ( I ) , I -J ,N  ) 

REALMS, 1009)  I  805  (  1 )  ,1  -1,  J  J ) - - 

J  *  JJ  ♦  l 
N  *  J  J  ♦  JJ 

REAOIS.  1000)  (  ...  BOBL  I  ) ,  1  -  J ,  N  ) - - — 

REAOIS, 1003)1  TPT2I I )  ,  1  >1,K  ) 

REAOIS. 1009)1  DPT ( I ) «  I  -1  <K ) 

00  S  I  *  l.JJ  - 

ETALGII)  -  ALPGIETAPTII)) 

N  *  JJ  ♦  I 

--  -  BOBIN)  *  ALOGIBOB(N)  )  - : - 

5  0051  I  I  *  ALOGIBOBI I ) ) 

REAOIS, 1009)  XN,  XZ»XEND»XPRINT»YO*DY 

-  N  *  XN  -  .  - 

NW  *  N/50 
1FINW.EQ.0)  NW  «  1 

-  ■-  NE2*  Nt>NE*. 00975  -  -  - - - 

VSR  =  OY  -  YO 
DY  -  VSM/XN 

51  =  EKPIYO) - -  - - * - 

NPl  *  N  ♦  l 
00  10  I  -  1 ,  N 

XI  -  FLOAT  (I  )  -  .5 - - - 

Y  *  YO XI /XN  *VSM 
EXPO  *  EXP(Y) 

PYEYII)  -  l./EXPO/OY  - -  - - - .... 

cTAI  *  c/PO  -  51 
50  *  t T  A l *♦ . 125 

2  T  A'Jo  ( I  )  x  BO*ETAI  ...  ..  -  .  -  . . 


NOT  REPRODUCIBLE 


IFIETAI  .CT. .00011  GO  TO  7 

Fi(l)  •  PPT ( l )  -  ...  - - - 

FK2)  ■  PPT ( 1 1  .♦  i) 

F 1 1 3  I  9  -.8*ET498m 

KF?4UU)  *  800(1)  ♦  B0B(JJ*1J  .....  _ 

.  GO  TO  a 

7  CALL  TEPP2N(II,3,ETAI  ,ETAPT,PPT,F1) 

XI  *  ALOGIETAI  l.  .  .  . . 

CALL  TERP2NIJJ.2.  XI  ,ETALG,BOB,Q) 

KF.VOIII  ■  0111  ♦  0(2) 

8  CALL  TcRP2N<  K»  1,.'1(2),TPT2«0PT,DPLS) - 

PH)  *  Fill) 

129(1)  ■  5809. /F 1I2)*DPLS 

.  E"P(  I  )  9  BO/P  ( 1 1  ...  - - 

DP(U$(I)  *  0PLS*2. 

OMlNUS(I)  -  234.*0PLU$( I ) 

ET*‘>3(IJ  9  ET  A98  ( I )  *COcL2*2.  _ 

T2TMII  ■  T29(T,*234. 

FP3Y(I)  *  Fl(3)/P(IJ*0YEY(i) 

-  -  XI  9  FI  ( 2 )  ♦  ♦.  5  - - 

KFNO(I)  ■  EXP(KFKO(  I )  -  32500./F1I2)  -  23.718998) /XI 
10  8(1)  »  NE2/Fl(2)/XI  /P(I)/P(I) 

00  20  I  9  1,51  - - 

xi  «  ;-i 

Y  *  YO  ♦  X 1/50.  *VSM 

EXPO  9  EXP(Y)  - - : - : - 

ETAO(I)  9  EXPO  -  Bl 
IFtI.eQ.l)  ETA3(1)  9  0. 

£TA90(I)  9  ET A9( I )**l. 125 _ 

IFd^'AQd  J.GT..OOOI)  GO  TO  27 
Fill)  9  PPT(1J 

--  -  FI (2)  ■  .PPT (  1 1  ♦  1) - 

r 1(2 1  •  -.8*ETa90( I ) 

GC  TO  28 

27  CALL  TcRP2N(  II,3,ETA0(I)  ,ETAPT,PPT  ,F1) _ .• _ 

28  CALL  TcRP2N(  K, 1,F1(2),TPT2,0PT,DPLS) 

PO(I)  9  Fill) 

...  FO(I)  9  FII3J  . - . . . 

FPO(I)9  Fl (3 l/Fl ( 1 ) 

DPLUO(I)  9  DPLS 

....  OPPIUd)  9  DPLS*234. _ 

ETA90I I )  9  ETA9Q(I)*C0EL2 
20  OTPLO(I)  9  H690./Fl(2)*DPLS 

- I  9  II  ♦  1  - - - 

XI  9  PPT(I)**.5 

90  9  Nc2/  PPT( I )/XI/PQ(  l  )/PQ(  l)  *4. 

...  KF.NOO  9  BCBUJ+l)  ♦  BOB(l)  .  _  _ _ 

KFNOO  9  EXP ( KFNOO  -  32500./  PPT ( 1 1  -  23.718998 ) /XI 
I  9  II  ♦  II 

.-.XI  9  PPT( I )**.5  .  ._  _  _ 

Bl  .»  licit  PPT(IJ/XI/PQ(51)/PQ(  51)  *4. 

CALL  TEKP2NI J J, 2 , 0. , ETALG, BOB, Q ) 

KFNUl  9  0(1)  ♦  0(2) 

KFNfJl  9  EXPtKFNOS  -  32509./  PPT(I)  -  23.718998J/XI 
RETURN 


IM1T* - EFN  .  SOURCE  STATEMENT — .-r — IFNIS)  — j*. 


SUBROUTINE  INITLZ  .  - - ...  -  - -  -  -- 

COMMON  ETAPTl  IDO  ),  PPT  (  300  l.ETALGUOO),  BOB  (200  ).TPT2  1100)  tDPTIlOO), 
1  0(2), FI (3) ,PHJ(401 ) ,NE2 , XN,  XZ • XENO, XPKINT , ENE ,0Y , YO, £T AI , XI , Y, 

?  =XPO,OPIS,A1,A2,A3,OELTAX,JOP, JOM,VSM,PHlY,DSYSP,DSYSM,SNH,  _ 

3  TSYSP.TSYSM 

CCw“3N/C0NSTl/XJ,QELJX, 2 J, 0Z25,DELJ 

C0MM0d/Cn.*JST2/Pl,P2,P3,P4,P5,Ps  . .  -- 

COXMON/CONST3/  DELZ ,  CDEt.2,  RHOECE,  VOX,  ENEE,  X212(KJ12 tOEL J12 ,0ELJX2 1 
1  UZ?  52 , C ( 8 ) , V fNEtCOEL, COELE, OOcCTA, OEtTA 

CO-XO.J/lNTLZ/P(IJ(401),SPJ(401),SMJtA01),JPJ(A01) ,JMJ{401),  . . 

1  USPJKCl  ),0SPJ!4jl  ),SPJH(400).SMJH(400),PHIJK(400) 

COMMON/  I.'IOl/  I, J ,M,  1 1 ,  J  J,K 

-  -  COMMON/ INOX/  IOX ( 6 )  - - 

COMMON/ NU.VJ/  N, NW , NP 1  ,fiC , NCOUN T 
REAL  Hit NE2 »JO°»JQM,JPJ,JMJ 

REAQIS, 10 )  PHIJ(1),P1,P2,P3,P4»P5,P6  _ 1 _ 

10  FORMAT  ( 7E10, 5 ) 

RcTOU,  11)  IDXI 1 ) «I 0X12) , 10X13), I0X(4)t 10X15), 10X16) 

_  11  F  C.l.VAT  1 615 )  _ _ 

XEND  *  XENC**.B 
CDEL2  »  C05L*C0EL 

t.'ic  «  l.6£-19*NE  .  .  _ _ _ 

ENc:  *  cNE/8  « 85E-14*,5 
CO;LE«  COEL/ENE 

DO  21  I  «  1,N  - - 

PHIJU)  «  PHIJ(I) 

PHlJHll )  *  PHIJ(l)  ♦  PHIJIl) 

-  .  srjiiu )  *  o..  ..  _ 

SMJHII)  «  0. 

SPJ(I)  *  0. 

- SMJ(I)  M.o. - : - 

JPJ(I)  «  0. 

J.MJ(I)  »  0. 

DSPJ(l)  »  0.  .  - - - . - 

21  OSMJI I  1*0.  .  • 

PHIJ(NPl)  *  PHIJIl) 

- SPJlNPl)-*  0. - : - 

SMJ(NPI)  *  0. 

JPJINPl)  *  0. 

.  JKJ1HP1)  *  0. - 

0S?J(NP1)  «  0. 

0SMJINP1 )  =  0. 

-  XJ  =  0.  - 

OEIJX  =  0. 

OELJ  *  0. 

0225  *  0. 

NC  =  0 

- N COUNT XZ/XPRINT - - - , - 

M  =  XZ 

OcLZ  *  XENO/XZ 

VOX  -  .2PV/0EU  - - - - 

OuElTa  *  COEL+OELZ 
RE  TURN 


NOT  REPRODUCIBLE 
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-STEP. - -- 


SOURCE  STATEMENT—^ _ IFNtS)  - 


SUBROUTINE  XSTEP  ....  - 

COMMON/CONSTl/X J.DEL JX, Z J«DZ25,0ELJ 

COMMON/ CONST 3/  DELZ « C0EL2.RH0ECE, VOX , ENEE, XZ 12, XJ 12.DELJ12.0El JX2 

1  0Z252  ,C(8 )  .V.NE.CDEL, CCEl.E,  DDELTA, DELTA.  _ _  _ 

COMMCN/INOl/  I,  J,M, II, JJ,K 
REAL  Ui 

XZ12  -  ZJ  ♦  DELZ  - - 

DELTA  *  CDEL*XZ 12 
DEL J 12  «  OELTA*OELTA 

-  XJ12  »  XZ12**1.20  .  - 

0ELJX2  *  DEL J12/XJ12 
OZ252  *  DELJ 12/ ( XZ 12**.25) 

C(  8)  »  DELJ 12  -  DELJ  -  - 

C(  l)  *  0ELJ12/RH0ECE 

C(  2)  *  VDX*(3.*OZ252  -  DZ25) 


IF(J.EQ.l)  dZ)  * 


C(  2I/.75 - 


C(  ?)  = 
C(  *)  * 
C(  S)  » 
C(  6)  » 
C(  7)  - 
C«  3)  » 
RETURN 
ENC 


V  *0=LJX2 
ENEE*05LJ12 

V  •  *(DELJX2  --DELJXI 
C(  8  l/RHOECE 

VOX* ( DZ2  52  ♦  DZ25) 
ENEE*C(  8) 


BOUND. - EFN  ...  SOURCE  STATEMENT — a — IFN I S 5 — — _ 


SU3R0UT 1NE  0CUNOO  - ...  - - - 

COMMON/ FUN3/  ET AQ ( 5 l ) , FQ( 51 ) , DTPLOI 5 l ) , FPQI 51 ) , PO ( 5 L I .DPLUOI 5 1) , 
l  0PM!Q(5l),ETA9Q(51) 

CGMl'.Q!'l/INTLZ/PHlJ(401),SPJ(40ll,SMJ(401),JPJ(40l)tJMJ(401)»  _ _ 

l  DSPJ(43l),DSMJ( 401  ),SPJH(400) , SMJHI 400 ) ,PH1 JHI400I 
COMfiON/MATl/X  (25  ),O(30?,  88(25) 

COM  •10.i7f-0NST3 /  OELZ,COEL2,RHOECE,VOX,ENEE,XZ12,XJ12,DELJ12,DELJX2 
l  OZ2»2,C( 8),V,NE,CDEL,C0ELE, ODELTA, DELTA 
COM.MON/CnNST2/Pl,P2,P3,P4,P5,Po 

COMMON/ CONST  I  /X  J  ,  DEL  JX ,  Z  J  ,  l)Z25  ,  DELJ  .  ...  .  _ _ 

REAL  JPJ.JMJ.NE 
X ( 3 )  =  10.E10 

X(  3)  *  10.E10  .  ..  - - - - 

0(25)  *  0. 

0(29)  »  0. 

RETURN  ....  .  ....  _ _ _ _ 

ENTRY  0011  NO  1 
X(  l)  »  10.E10 

X(  7)  a  10.E1C - 

X(2b)  -  10.E10 
0 (2o )  a  0. 

0(27)  a  0.  - 

0(30)  a  P1*OOELTA*X ( l ) 

RETURN 

-  E  NO  - 


EMPTY...—  - .  _CFi4  SOURCE  STATEMENT _ - — IFNIS)—- - 


SUBROUTINE  EKPTY  - - - - - - — 

COY-O.N/CONST3/  DELZ  ,C0EL2*RH0ECE ,  VDX,ENEE,X£12,XJ  12,0E1.J12  »DEt_  JX2 , 
1  DZ252,C(8) ,V,N£,CDE(.,CDEIE, ODELTA,  DE1.TA 
.  CCPKON/FUNI/CPIUSKDO)  ,OM1MUS(40D>  ,EPI400) ,ETA98I400) .T29(400)»  - 
1  B(40jJ  .FP'YUDO  J,KFNO(>.ODJ,DYeYU00  5,P(4CO),T29KC400) 
COMKO.J/1NTLZ/PH1  J  (401),  SPJ(  4011,  SMJt4DH,  JPJ(  401  »,JNJ  1401)  , 

.1  USL,J('tOl),DSHJ(4Gl),SPjH(400),SHJH(  4Du )  ,PKIJH(40C) . 

COwKUN/HATl/X (25 ), 0(30), 00(25) 

COMMON/ I N01/  1, J,M, I !,JJ,K 

.  CD-MO-i/INOK/  10X  ( 6 )  -  -- - - - 

REAL  Nc,  JPJ,  JHJ.KF.NO 

OD  sG  11  »  1,25 

-  BO*  11)  *  0.  ..  - 

x ( 1 1 )  «  0. 

50  0(11)  »  0. 

RETURN  ...  -  - - - — 

ENTRY  FIRST 
K  «  1-lOXIl) 

D(2a)  «  00ELTA*ETA98(K) .  . . . . 

X(ll)  «  DPLUSIO  ♦  0ELTA*ETA98(K) 

0  (  2j  )  «  0(26 )/( X ( 11)  -  0(26)) 

..  X(l)  «  T29IK  )/XI  11)  - - - 

X  ( 21 )  *  X ( l ) ♦ (SPJHI K )  ♦  OSPJIK+l)  ♦  OSPJ(K)) 

X(l)  «  X(ll*PHlJH(Kl 

-  X « 1 1 )  «  -P(K)/X(ll)  - - 

0(2t>)  «  I  XI 1 )  ♦SPJHI  K)  ♦  XI 11  )*(  JPJCK+i)  *JPJ(  K) )  1*0(26) 
IFII'JXIIKEQ.I)  SO  TO  80 

...  D(l)  «  X(l)  ♦  DYEY(K)  . . . . . . 

XIII  «  XCl)  -  OYeY(K) 

0(11)  «  X  1 11  ) 

0(21)  -  X 1 2 1 )  - 

RETURN 

80  00(1)  »  X(!)  -  OYEY(K) 

..  .  XII)  *  X(l)  ♦  DYEY(K)  — - 

88(11)  >  X(ll) 

06(21)  «  X 1 2 1 ) 

-  RETURN  -  - - - 

ENTRY  SECOND 

K  *  I  -  10X12) 

_ 0(27)  •  D0ELTA*ETA98(K)  _ 

X 1 1 7 )  •  DMINUS(K)  ♦  0ELTA*ETA981K) 

0(27)  «  0(27 I/IXI17) ”01 27)) 

_ X I  7 )  »  “T 29M ( K ) / X ( 17)  - - 

X ( 22 )  *  XI7 ) ♦ (SPJHI K)  ♦  OSMJIIUI)  ♦  DSMJ(K)) 

X  (  7 )  «  XI 7 ) *PHl JHI K ) 

- X ( 1 7  )  «  -PIKJ/XI17)  - - - : 

0 1 27 )  «  IX(?)*SHJM|K)  ♦  XI 17 ) *  I JMJ(K  +  1)*JMJIK)))*0(27) 

I F  t IDXt  2) #EQ • 1 )  SO  TO  10 

-  0(7)  a  X ( 7 )  ♦  OYEY(K)  - 

X 1 7 )  *  X  1 7 )  -  DYEYIKI 
0(17)  »  XI 17  ) 

0122)  «  X(  22  )  ..  .  - - - 

RETURN 

10  B H ( 7 )  a  X ( 7 )  -  OYEY(K) 

X(  /)  a  X 1 7 )  ♦  DYEYIK)  .  .  ...  . -  - 


36 


88(17)  ■  X( 17 ) 

--  83(22)  a  X(22) -  -  -  ...  .  - - - - 

RETURN 
ENTRY  THIP.D 

K  *  I  -  I0X(3)  .  .  ..  .  - - - 1 

XU)  *  3(K)*C(1) 

XO)  3  -X(8)*SMJH(K)  -  C(2)*EP(X) 

X(i)  *  -X(9)  *  (SPJH(K)  ♦  DSPJU  +  1)  ♦  DSPJ(X) ) - 

X ( 13 )  *  C(3)*F?EY(K) 

D(2i)  =  C(5)*FP£Y(K)*!SPJ(K»1)  -  SPJ(K))  *C ( 6 > * (-KFNO( K)  ♦  8(K)* 

1  SRJH(K)*SMJH(K) )  -C(7)*EPm*(OSPJOX*l)  ♦OSPJ(K))-  -  - 

I F ( I OX ( 3 > .EQ.l)  GO  TO  20 
■0(  Ji  3  X  ( 3 )  -  X  (  13 ) 

X  (  3)  *  X ( 3 )  ♦  X  ( 13 )  - - 

0(<S)  *  X  (  8  ) 

0(13)  =  OYEY(K) 

X(lJ)  =-OYEY  (X )  .  .  ......  _  - - - - 

RETURN 

20  Btt(j)  »  X ( 3 )  ♦  X ( 1 3 ) 

-  X  ( 3 )  =  X(3)  -  X  ( 13)  ... _ 

88(3)  *  X(0) 

83(13)  3-0YEY(K) 

X ( 1 3 )  *  DYEY(K)  _ _ _ _ _ 

RETURN 

ENTRY  FU'JRTH 

-  k  -  I  -  i nx( 4 1  -  - _ - 

X (4 )  3  B(  K ) *C ( 1 ) 

X(9)  »  -X(4)*SPJH(K)  -  C( 2 1 *EP( K > 

.  X  ( 4 )  *  -X(4)*(SMJH(K)  ♦  DSMJ  ( <  ♦  1  >  ♦  OSHJ(Kj) -  ....  _ 

X  ( 1  7 )  *  C  (3 )  * FPE Y  ( K ) 

D(2v)  *  C(5)*FPEY(K)8(SHJ(K*1)  -  SHJ(K))  +C ( 6 ) * ( -KFNO( K)  ♦  B(K)P 

...  I  SPJH(x)«SMJH(K)  )  -C(7)*cP(K)*iDSMJ(KU)  DSXJIK.)  ) - 

IF(IDX(4) .=0.1)  GO  TO  30 

U(4)  *  X ( 4 ) 

-  0(9)  *  X (9 )  -  X(19>  . . .  ...  ..  ... - 

X ( 9 )  *  X l 9)  ♦  X ( 19 ) 

0(19)  »  DYEY  ('< ) 

- X  ( 19 )  s-DYEY  (K  )  ...  .  - - - - 

RETURN 

30  BB(4)  •  X (4 1 

.  .  03(9  )  *  X  (9  )  ♦  X ( 19  )  - - - - - - 

X  ( 9  )  *  X  !  9  )  -  X  (  19  ) 

BB ( 1 9 )  3-OYEYU) 

-  X  ( 1 V )  =  DYEY(K)  _.  - - - 

RETURN 
ENTRY  FIFTH 

r.(5)  =  C(4»/p(x) 

X(10)3  -X (5 ) 

0(33)  =  -C(8)/P(K)*(SPJH(K)“SMJH(K) ) - 

I F ( I  OX ( 5 ) . EQ . 1 )  GO  TO  40 
0(5)  =  X ( 5 ) 

0!  10)  *  x(10  )  . . . . .  . 

0(2a)  *  DYEY(K) 

X(25)  3-OYEYIK) 

RETURN  ....  ..  ...  .  ...... 

40  88(5)  *  X(5) 

88(10)  *  X(10) . . 

*0(25)  3-DYEY«) 

X (25 )  =  OYEY(K) 

100  RETURN.  .  ........  _ 

END 
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K£ST *  - -  -  EFN  SOURCE  STATEMENT.—-- — IFN1S) — -- 


SUBROUTINE  REST (Z I  -  -  - -  - 

DIMENSION  Z1  103.00) 

COMMON  £T APT ( 130) ,PPT1 300) .ETALGI 100) .038(200) .TPT21 100) tOPTI 1001 , 

1  Q(2).F1{3),PHJ(401).NE2.XN,XZ,XEND,XPRINT,£NE,0Y,Y0,ETAI,XI,V,  - 

2  E<PO,OPLS.Al,A2fA3,OELTAX,JOP,JOM,VSM,PHIY,DSYSP»DSYSH,SN.H, 

2  TSrSP.TSYSH 

COm.U/COUST l/XJ,DELJX,ZJ,0Z25,DELJ  .  . 

CU*O'.0.4/C0NST2/Pl,P2,P3,P4,P5,P5 

CO^/.U.'4/COUST  3/  DELZ  , CDEL2 , RHOECE, VOX ,ENEE , XZ 12?X J12 »OEL J12 ,OEL JX2 » 
1  OZ2S2, CO). V,SE,COEL,CO£L£,OOELTA, DELTA  -  .  -  .  -  - 

CUw;’.0h/FUN1/P  PLUS  1400) , DMINUS1  400  )  » EP14CD)  •  5TA98C  AOO)  .T29C400 )» 

1  0(<.0u)  ,FPEY(  430  ).KFHO(  400),  DVEYC  400 ),  PC  400).  T29M  (400) 

-  euv»0:i/ru.‘l3/  EfAO(01)iFO(al)tOTPLO(51)iFPO(51)tPQ(51) .  DPLUOt  51)  »  — 
1  UPHlQ(Sl)tETA3C(3l) 

COXRC.R/IMTlZ/PHI  J1 401).  SPJl  401),  SMJ1 401  J.JPJ  1401),  JMJ1401)  , 

1  DSPJ ( 43 1 ) , DSMJ ( 40 1 j , SPJH1 400),  SMJH1  40u) .PHI JH 1400) -  - - 

C0mw/HAT2/XX120D:>) 

CU»KK(/IN01/  I.J.M.I I, JJ.K 

.  -  COF’fiu.'l/ 1 'IDX/  10X16)  -  - - - 

CO* '.On/ MUHB/  N, NW ,  NP l ,  UC , NCOUNT 
REAL  JOP.JOH.KFMO, JPJ, JKJ.NE ,NE2 

.  READY  FOR  BACK  PJLTIPLICAT  ION  .  _ _ , - - 

JJ  » 

1 1  *  23*N 

--  00  2iO  1  ■  1,N  . . . 

11  »  11  -  25 
JJ5  *  JJ 

jj  »  jj-5  - - 

00  250  X  *  1,5 
1.1  *  JJ  ♦  K 

...  -  1.2  *  11  ♦  K - 

00  250  L  *  1,5 
L3  »  JJ5  ♦  L 

<X( L 1 )  *  XX(Ll).-  Z1L2)*XX1L3>  - 

250  L2  «  LZ  ♦  5 

ANSWs»S  ARE  IN  XX  ARRAY.  WRITE  ANO  UPDATE  BEFORE  STEPPING  X. 

-  -  K  *  -4  - . -  —  - - 

UO  2oO  1  -  l.NPl 
K  *  K  «-5 

.  .  .  OSPJ(l)  *  XX(K)  .  _ - -  - - „ 

OSKJII)  *  XX 1 K* 1  ) 

spjm  *  spj(i)  ♦  ospjii) 

_ SMJ(I)  «  SMJ1I)  ♦  OSY.JII) - l — 

JP  )( I )  *  JPJ(l)  ♦  XX ( K+2 ! 

JMJ( I )  -  JMJ(  I)  ♦  X  X 1  K 1 3  ) 

PHIJ(I)  *  PH  I J I  I  )  ♦  XX(K»4)  _  —  .... - 

IF1I.E0.1)  CC  TO  260 
M  =  1-1 

S P JH 1(1)-*  SPJ1I)  ♦  SPJ1II)-- - 

SHJH(II)  =  SPJ1I)  ♦  SMJ(Il) 

Pill  Jill  (1)  *  PH!  J  1  I )  ♦  PHIJl  III 

260  CONTINUE  - - : - - — 

XJ  *  X  J  12 
UELJ  =  DcL J l 2 

UELJX  =  0ELJX2  . •- -  - - 


NOT  REPRODUCIBLE 


I J  »  XZ 12 

-0Z25  «  3,2152 - - - 

NC  *  NC  ♦  I 

IFI’iC.LT.NCOUN? )  RETURN 

C  WRITE  HEADING 

WAITcIS.lDDl)  XJ.OELZ.N 

1001  FORMAT (5H1X  «  , F7.2, 5X, 10HDEITA . Z  .»  .F9.4.SX. 16KETA- INTERVALS  «  * _ 

1 1 3  */8dX  .5HKP LUS  » 4X»  6HU  PLUS.9H  0  H I NUS « 32. 1 6H0TPLUS , 9H  OTMNUS  / 

23  X,  _;Hi  T  A*  5X*  6HN.M  INUS  »4X  *  5I<SPLUS»4X»  6HSM INUS.  5X  ,  5HJPLUS  ,  4X,  6H  JR  I  N'US 

3  >4X,oUPHIETA,5X,3HPHI,7X,2H-V,  5X,5H*PH1V,4X,;h*jY/S.4X,5H*SY/S.  _ 

44 X  .5M*SY/S  *4X.5H*SY/S  //$ 

C  COKPUTf  PHI  I  «  PHJ) 

-  IFIIUXIaI.EO.il  CO  TO  231  -  —  -  - —  .  — - - - 

PIIJI.4P1)  »  0. 

00  232  1  ■  l.N 

-  _  K  *  NP1  -  1  _  _ _ _ 

252  PHJIK)  »  PHJiK+lt  -  .54PHI JHI X ) /OYEYIK I 
GO  TO  275 

-  251  PIW(l»  ■  0.  -  -  -  ■  -  -  - - - - 

DO  253  1  ■  l.N 

253  PlIJII  +  l)  «  PHJ(I)  ♦  ,5*PHl JHC  il/DYEYCll 

...  275  A I  «  V*DEL  JX  .  .  - - - 

A2  »COELE*ZJ 
DELTAX  «  -DELTA/XJ*V 

11*0  - -  - 

DO  :0G  l  •  1.NP1.NW 
11  *  11  ♦  1 

.  aj  «  AiPFPonn  - - - 

JUP  «(JPJ(I)  -  A3*SPJ(1»)/A2 
JOY  ■ ( JHJ II)  -  A3*SHJ ( I >  J/A2 

-  A3  *  PUI 1 1 ) /DELTA  .  .  - : - - 

VSP  «  0ELTAXYF3I III 
SUM  «  Ne*$M  HII/POIII) 

.  PHIV  a  OTPLOMI  I *PHI J (  1 1/OELTA  .  - - - 

DSYSP  »  “PHI Y  ♦  A3*JPJ( Ii/SPJI I J 
USVSH  •  PHIY*234.  ♦  A3*JMjm/SHJIII 

- A3  *  ETA90II I l*OELT A  - - 

TSYSP  »  OPLUCIIIl  ♦  A3 
TSVSH  *  DPMI C ( I  I  1  ♦  A3 

—  .  DSYSP  ». 0SYSP*0PLU0( 1 1 I/TSYSP - - - 

DSYSM  *  OSYSP*nPHIQi; II/TSYSM 
TSYSP  *  0SYSP*A3/UPLUQI III 

- TSVSH  a  DSYSK*A3/0P«I0< III  - - 

IFISPJIII.NE.O. »  GO  TO  280 
DSYSP  *  1.E38 

-  TSYSP  *  I.E3R  _  _ 

280  IFlSv.jm.NE.O.)  GO  TO  300 
DSYSM  »  1.E28 

T SYjM  s  1  .S3B  -  - -  - 

303  WRITE  (6*.  1 333  I  ETAC(II).  SNM.SP JI  II.  SHJ  (  I )  . JOP. JOM. PHI  Jill  .PKJI 1 1  > 

1  VSK, PHI Y, DSYSP. OSYSM, TSYSP. TSYSH 

-  1003  FORMAT 11P2E9.2.1P5E10.3.1P7E9.2)  _  _ _ _  _ _ 

IFIXJ.GT.P51  STOP 

RETURN 

ft  *i n  * 

.  -  -  - 
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.  XOCOR 


. - EFN  .  SOURCE-  STATEMENT — —  IFN(S) 


SU3*OUTlNE  XOCORR(Z)  -  -  —  - 

OlVEMilUN  Z ( 10300  ) 

C0HMGN/HAT1/X (25), 0(301,00(25) 

CGVv.0N/KAT2/XX(  200S )  - - - 

CCV.-O.-l/lN0l/  I.J.K,  IltJJtK 
COV-ON/1NOX/  IOX  (  6 ) 

JJ  *  25*(  1-2)  - - 

JJ3»  3*1  -  9 
JJl  *  JJO*  l 

JJ2  *  JJO*  2  .  .  .. - - - - - - 

JJJ  *=  JJO*  3 
JJ-.  *  JJO*  4 

—  I  F  t  IUX(  D.EO.O)  CO  TO  7 - 

00  3  II  «  1.21.5 

K  ■*  JJ  ♦  II 

5  XIIll  «  X(!I)-BB(1)*Z(K  >-BB(  ll)*Z(K*2)-RB(21)*ZIIt*4)  - 

0(2v)^0(26)-3BI 1)*XXI JJ0)-BB( 11 ) *XX( JJ2 )-BB ( 21 )*XX( JJV) 

7  IF(IUX(2J.EQ.0>  GO  TO  12 

00  lw  II  «  2,22,5  - - 

K  *  JJ  ♦  II 

10  XIII)  »  XIII  )-9B(  7)  *Z(K  )-BBU7)*ZIK*2l-BB(22)*ZIK*3) 

0(Z7)*0(27)-BB(  7J*XXUJ1)-BQ117)*XX(  JJ3 )-BB( 22)*XXI  JJ4) - 

12  IF(IOxm.EO.O)  CO  TO  17 
00  13  II  <•  3,23,0 

X  :  JJ  HI  . . . . 

15  XIII)  «  Kill  )-BBm*m-2)-BB(  8 )*Z( K-l  )-B8(  13)*Zf X) 
OI23)=OI23)-BB(?)*XX( JJO)-BB<  B )*XX« JJ 1  )-B»( 13)*XX( JJ2) 

17  I F  {  I'JX(4l.E0.C)  GO  TO  22  -  ....  .  . . 

00  2 £  II  «  4,2*. 5 
K  «  JJ  ♦  II 

20  X  (  1 1  )  *  XUI)-BBI4)*Z(K-3)-BB(  9)*Z(K-2)-BBC  19)*Z(K)^ - 

0(2‘;)=0(29)-P.B!4)*XX(JJO)-BBI  9)*XXl JJl  )-BB I 19 ) *XX( JJ3) 

22  IF( 10X15). EO.O)  RETURN 

00  25  II  «  5.25,5  -  .  . . . 

K  *  JJ  ♦  II 

25  X ( 1 1 1  *  XIII )-BB(5)*Z(K-4)-BB(10)*2(K-3)-Bff(25)*7tK) 

- D(5C)*0(30)^BB(5)*XX(JJ0) ~BB( 10)*XX( JJl )-BBt  25)*XXIJJ4) - 

RETURN 

END 


. .  ENOO. - :_-_EFN_  SOURCE  ST2T EMENT- — - — JFNCS) — t _ 


subroutine  enod  . . 

C0W0N/1N0X/  .10X16) 

IF(IOXIl).EQ.O)  CALL  FIRST 

I F ( 10X12)  .EO.O)  CALL  SEC0N0 - 

I F 1 1  OX  t 3) .EO.O)  CALL  THIRD 
IF ( 10X14) .EO.O)  CALL  FOURTH 

IF(IDXC5).EQ.O)  CALL  FIFTH  -.  - 

CALL  BUUNOO 
RETURN 

ENTRY  ENOl  _ 

IF(IOXm.EO.l)  CALL  FIRST 
I F( IOXI 2) .EQ. 1)  CALL  SECOND 

IF(  10X13). 50.1)  CALL  THIRD _ 

I F ( IUX(4) .EQ. 1)  CALI  FOURTH 
IF ( IOXI 5) .EQ. 1)  CALL  FIFTH 

CALL  BOUNOl  _ 

RETURN 

END 
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